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, ABSTRACT 

Adaptive optics observations of the flattened nuclear star cluster in the nearby 
Q . edge-on spiral galaxy NGC 4244 using the Gemini Near-Infrared Integral Field Spec- 

■ trograph (NIFS) have revealed clear rotation. Using these kinematics plus 2MASS 
I photometry we construct a series of axisymmetric two-component particle dynamical 
. models with our improved version of NMAGIC, a flexible x^-made-to-measure code. 

The models consist of a nuclear cluster disc embedded within a spheroidal particle 
population. We find a mass for the nuclear star cluster of M = I.GIq ^ x 1O''M0 within 
^ ■ ~ 42.4 pc (2"). We also explore the presence of an intermediate mass black hole and 

show that models with a black hole as massive as M, = 5.0 x 10^ Mq are consistent 
\^ . with the available data. Regardless of whether a black hole is present or not, the nu- 

clear cluster is vertically anisotropic {f3z < 0), as was found with earlier two-integral 
\ models. We then use the models as initial conditions for iV-body simulations. These 

QQ ■ simulations show that the nuclear star cluster is stable against non- axisymmetric per- 

■ turbations. We also explore the effect of the nuclear cluster accreting star clusters 
Ps| , at various inclinations. Accretion of a star cluster with mass 13% that of the nu- 
clear cluster is already enough to destroy the vertical anisotropy, regardless of orbital 
inclination. 
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1 INTRODUCTION 

Studies of the centres of galaxies across the Hubble sequence 
have shown that they frequently host central massive objects 
such as massive nuclear star clusters (NCs) and supermas- 
sive black holes (SMBHs). NCs are present in roughly 75% of 
low and intermediate luminosity dis c and elliptical galaxies 
l|Boker et al.ll2002l : [Cdte et al.ll2006l '). These NCs are intrin- 
sically very luminous, with typical Mj ~ —12, and size s 
similar to globular clusters (re// ~ 5pc: lBoker et al.lboOll '). 

Two hypotheses have been offered to explain NC for- 
mation. One scenario envis ages NCs forming in situ out of 
gas cooling onto t he centre (|Milosavlievidl20o3 : iBekki et al.l 
l2006l : iBekkil 120071 ). Alternatively, NCs may form from star 
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clust e rs merging at the centres of galaxies (Trcmain e et all 
1 19751 : iLotz et all [20011 : ICapuzzo-Dolcetta fc Miocchll2008l '). 
Which hypothesis is correct determines whether NC growth 
is limited by the supply of star clusters from the host 
galaxy or regulated by fe edback from in situ star formation 
(McLaughlin et aLllioO^ V 

NCs exhibit several scaling relations. The l uminosity of 
NCs correlates with that of their host galaxy (|Boker et al.l 
l2002l : IC6te elaLlbOOdlErwin fc Gadottfeoiol v A number of 
studies also found that their mass, Mnc , correlates with the 
velocity dispersion of t he host bulge, the Mnc — a relation 
(iFerrarese et al.l l2006l : IWehner fc Harris! l2006l : iRossa et all 
I2OO6I ). Early work found that this Mnc — c relation is par- 
allel to the M. - a relation of SMBHs (jCebhardt et all 
l200d IFerrarese fc MerrittlbOOOl ). with NCs being about 10 x 
more massive, at the same a, as SMBHs. However recent 
work has questioned how comparable NCs and SMBHs are. 
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lErwin fc Gadottil l|2012ah find that NCs and SMBHs fol- 
low different relations, with SMBH masses correlated with 
properties of the bulge, while NCs seem to correlate bet- 
ter with propertie s of t he entire host gal a xy. In stead, both 
iLeigh et al.1 and IScott fc GrahamI l|2012l ) show that 

there is an Mnc — o" relation but with a significantly differ- 
ent slope than for SMBHs. Part of these differences may be 
due to the fact that t he M. — a and Mnc — c relations de- 
pend o n Hubble type (jCreene et al.ll2010l :lE rwin fc Gadottil 
l2012bl '). and especially that SMBHs and NCs in late- type 
galaxies not following the same scaling relations as in early- 
types. It is therefoere important to shed light on how mass 
accretes into the centres of the late- type galaxies. 

It is generally thought th at AGN f eedback is responsible 
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narios where this rel ation arises because the galaxy regulates 
SMBH growth fe.g. [Surkert fc sii^l200ll: iKazantzidis et all 



I2OO5I : iMiralda-Escude fc KoUmeieij 120051) or purely indi 
rectly |Hachnclt fc Kauffm annl I2OOOI : lAdams et all l200ll . 
I2OO3I : [jahnke fc M accio 20l3) have also been proposed. 
Therefore NC scaling relations hold the promise of illumi- 
nating the origin of SMBH scaling relations. The similarity 
of the Mnc — o" and M, — a relations suggests that related 
processes may be at wo rk in controlling their growth (e.g. 
iMcLaughlin et al.ll200^ . 

So m e galaxies host both a NC and a SMBH (|Seth et al.l 
l2008al lbl: [Graham fc Spitleill2009l ') and their properties pro- 
vide further constraints on the relation between the two. 
For instance iNeumaver fc Walcheil (|2012l ') found that NCs 
and SMBHs do not coexist at M./Mnc > 100, which they 
propose suggests that NCs are precursors of SMBHs. Con- 
structing an (M,-|-Mnc) — f" relation which includes the mass 
of both the N C and of the SMBH (Graham ct al. 2011), 
iGrahamI ||2012D found a flatter relation than the M, — a re- 
lation. 

Unlike SMBHs, the assembly history of NCs can be 
constrained from their morphology, stellar populations and 
kinematics. In late- type spirals, NCs have been found to con- 
sist of multiple stellar populations, typically a young popu- 
lation ( < 100 Myr), and a dominant population older than 
1 Gvr (iDavidge fc Courteaul |2002|: ISchinnerer et al.1 12OO3I : 



iRossa et al.1 12006*: 'Wa lcher et all 120061 ). The Huhble Space 
Telescope has revealed that the NCs of edge-on galaxies host 
multiple stellar popula tions associated with different mor- 
phological components (ISeth et al ]|2006l ). These NCs consist 
of young blue nuclear cluster discs (NCD) and older nuclear 
cluster spheroids (NCS). Optical spectra of the NC in the 
edge-on Scd gal axy NGC 4 2 44 (i « 90° ) , the nea rest galaxy 
in the sample of lSeth et all (|2006D (D=4.37 Mpc: ISeth eratl 
l2005ah . indicates the presence of multiple stellar popula- 
tions, while near in frared spectroscop y showed that the NC 
is rapi dly rotating (ISeth et al l2008bl ). Using A'^-body simu- 



lations |Hartman^^et^|al] 1 2011 . hereafter Hll) showed that 
the NC in NGC 4244 cannot have assembled more than half 
its mass via the accretion of star clusters. Hll also mod- 
elled the NC in NGC 4244 using two-integral JAM models 
l|Cappellarill2008l ). obtaining a mass of (1.1 ±0.2) x lO'^ Mq. 

To help shed light on the formation of NCs in late-type 
galaxies, in this paper we study the NC in the nearby Sc 
galaxy NGC 4244. The outline of this paper is as follows. 



Section[2]describes the observational data and how they are 
used in the dynamical modelling. Our modelling method , 
the x^-M2M code nmagic based on lde Lorenzi et al.) (|2007| . 
i2008, . 2009), is described in Section |3] including additional 
code development. We construct various axisymmetric par- 
ticle models of the nuclear region of NGC 4244 in Section |3] 
using this improved code. The models consist of a NCD and 
NCS having separate mass-to-light (M/L) ratios. Using the 
best model as initial conditions for A''-body simulations, we 
explore the evolution of the NC in Section [5] Section [6] dis- 
cusses our results in the context of NC formation. 



2 OBSERVATIONAL DATA 

We begin by describing the photometry and how these data 
are deprojected to obtain a three dimensional luminosity 
density. After this, the integral-field kinematic data are pre- 
sented. We adopt a distance to NGC 4244 of 4.37 Mpc 
ISeth et al.l l2005bl ). At this assumed distance, 1" corre- 
sponds to 21 pc. 



2.1 Photometry 

Here we give a brief summary of the photometric data and its 
model representation, both described in detail in lSeth et al.l 
(2005a, 2008b). 

The photometry consists of K-h&nA data either from 
2MASS or from the NIFS observations of ISeth et al.1 
(2008b). The main disc (MD) of NGC 4244 hosts at its cen- 
tre a NC which is composed of a NCD and an oblate NCS. 
The Jf-band mass-to-light ratio M/ Lk is estimated to be 
0.5 — 0.75 for the galaxy as a whole (from integrated colours 
take n from LEDA \B - V = GA- 0.6] combined with M/L 
from iBell et"al] (|2003t )). T he M/Lk of the NCD is in the 
range 0.1 — 0.25 based on iBruzual fc C hariot (2003) mod- 
els applied to the optical spectroscopy of Seth et al,, (200(j ). 
which agrees well with the fitted luminosities for the disc 
in HST/ ACS and NIFS bands. The NCS stellar populations 
are poorly constrained and a M/ Lk between 0.5 and 1.2 is 
likely. 

The surface brightness of NGC 4244 is decomposed 
into an axisymmetric three component model. The lu- 
minosity distributions of the MD and the NCD are 
modelled as projected edp e-on exponential discs {e.g. 
Ivan der Kruit fc Searlj|l981 



T.{x,z) = Eo ( ^ ) Ki{x/hr)sech{z/zQ)^ 



(1) 



where Eq and hr are constants, and K\ is the modified Bessel 
function. The corresponding model parameter s for th e MD 
and NCD models are taken from ISeth et al.l (|2005al ) and 
ISeth et al.l (|2008bh . respectively. 

On the other hand the NCS is represented using a lSersid 
(,196&) profile: 

I{x, z) = h exp(-6„((i?/i?e)'''" - 1)), (2) 

where Je is the surface brightness at the effective radius 
Re, bn 1.992n - 0.3271 and R = y/sF+JT/q)^ is the 
elliptical radius, wit h flattening q. Bes t fit parameters have 
been obtained from lSeth et al.1 l|2008bl ). 

All models were convolved with a Gaussian point spread 
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So 
[L0/pc2] 


PO 

[Lo/pc^] 


hr 

[pc] 


zo q 

[pc] 


[L0/pc2] 


Re 

[pc] 
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MD 


598 


0.167 


1783 


469 








NCD 


1.41 X 10^ 


2.08 X 10* 


3.39 


1.19 








NCS 








0.81 


8.73 X 103 


10.86 


1.68 



Table 1. Best-fit parameters of the photometric three component model taken from lSeth et al.l l l2005al . l2008bl V The model consists of 
an exponential main disc (MD), an exponential nuclear cluster disc (NCD) and a Sersic nuclear cluster spheroid (NCS). 



function (PSF) of 0V23 FWHM during the fitting process, cf. 
ISeth et al.l(|2008bl ). The best fit parameters are summarised 
in Table [T] The left panel of Figure [T] which presents the 
NCS model, shows its surface brightness, hk, and ellipticity, 
e = 1 — q, profiles. 



2.1.1 Deprojection 

To compute the three dimensional luminosity distribution, 
each component of the surface brightness model is depro- 
jected individually. The edg e-on deproject ion of an ax- 
isymmetric system is unique (|Rvbicki|[T987l '). The surface- 
brightness profile of Equation [T] corresponds to an expo- 
nential disc, so the deproje ction is readily given by {e.g. 
Ivan der Kruit fc Searlelll98ll ): 

p{R, z) = po exp(— sec\i{z/zof, (3) 

where po = So/(2/!,r), and Eq, hr and zo are as in Equa- 
tion [T] Their values are given in Table [T] 

Unlike the surface brightness profile of an exponential 
disc, the Sersic profile of the NCS cannot be deprojected 
in clos ed form. We therefore use the program of iMagorrianI 
l|l999[ ) to numerically deproject the surface brightness distri- 
bution of the NCS. The program finds a smooth axisymmet- 
ric density distribution consistent with the surface bright- 
ness distribution for the specified inclination angle (here 
edge-on, i.e. 90°), by imposing that the solution maximises a 
penalised likelihood. Because the deprojection is computed 
numerically and tabulated on a grid, the reprojected surface 
brightness profile may not match the Sersic one perfectly. A 
comparison of the NCS Sersic photometric model and its 
edge-on deprojection reprojected onto the sky plane, seen 
in the left panels of Figure [1] shows that the numerical de- 
projection is in fact very reliable. The right panel presents 
iso-density contours in the meridional plane of the NCS ob- 
tained with ,Magorriani (,199a )'s code. 



2.2 Kinematic Data 

The integral-field NIFS kinematic data we use in the mod- 
elling were presented in ISeth et al] (|2008bl ) and consist of 
velocity, velocity dispersion and the higher order Gauss- 
Hermite moments /iS and /i4 (Ivan der Marel fc Franx|[l993 : 
lGerhard|[l993 '). The NIFS field-of-view extends to ±1.5" 
along each direction, but the usable data is within ±1.0". In 
this field-of-view, the positions of the "spaxels" within which 
spectra were taken define a grid of 63 x 71 = 4473 cells, 
which serves as the basis grid for Voronoi bins for which the 
velocity, velocity dispersions, /i3 and /i4 are given. 




Figure 1. Left: Comparison of the NCS Sersic photometry with 
the reprojected three dimensional luminosity model. The data 
points correspond to the NCS Sersic profile, which was derived 
via a morphological fit. The solid line shows the model projected 
onto the sky plane. The upper panel shows the surface brightness, 
11 and the lower one shows the ellipticity, e = 1 — q. Right: 
Iso-density contours of the NCS in the meridional plane. The 
luminosity density was obtained by an edge-on deprojection of 
the NCS Sersic surface brightness. The contour labels are given 
in units of logjQ Lq pc~^. 

3 METHODS 

We construct a range of dynamical models for the NC of 
NGC 4244. These models consist of a disc and a spheroidal 
particle population representing the NCD and NCS re- 
spectively. We use an adapted version of the flexible x^- 
made-to-measure (M2M) par ticle code nmagic described in 
Ide Lorenzi et al.l (2007I . l2008h . This section describes a few 
ingredients required to construct dynamical M2M models 
and presents further development of NMAGIC comp ared with 
our previous work in Ide Lorenzi et al.l (|2007l . 120081 1. 

3.1 Model observables 

The central luminosity volume density of the N CS is about 
a milli on times larger than that of the main disc jSeth et al.l 
l2008bl ). Thus the NCS dominates the luminosity distribution 
in the central region out to the edge of the observational 
data, allowing us to neglect the luminosity distribution of 
the main disc when computing the photometric observables. 
At a distance 7? = 30 pc from the centre, the luminosity 
density in the equatorial plane of the NCS is still about 100 
times that of the MD. 

We compute separate spherical harmonic coefficients 
Aim for the density of the NCD and the NCS. The cor- 
responding errors ar e inferred followin g a Monte-Carlo pro- 
cedure described in Ide Lorenzi et all (2008), in which the 
AimS axe computed many times from random rotations 
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about arbitrary axes of a suitable particle realisation. Here, 
we used the isotropic particle model generated from the 
major-axis density profiles of the NCS and NCD compo- 



nents, described in Section [3.21 below. The A 



and Ai 



constrain the photometry of the disc and spheroidal par- 
ticle populations, respectively. We use even Aim^s up to 
Imax = 8 in 40 radial bins, unevenly spaced, for a total 
of NAim = 2000 photometric constraints. The grid starts 
at r™i„ = O'.'OOOl (0.002 pc) and extends to r^a^ = 2" 
(42.4 pc). 

The NIFS kinematic data (velocity v, dispersion a, h3 
and hi) are bi-symmetrised by adopting a point-symmetric 
reflection with respect to the cen tre of the galaxy, as de- 
scribed in Ide Lorenzi et al.l (|2009l '). followed by a reflection 
about the major axis. The resulting data within ±0.7" are 
shown in the upper panel of Figure [SI As kinematic ob- 
servables we use lu minosity-weight e d Gauss-Hermite coeffi- 
cient s hi up to /i4 ( GerhardI 1 19931 : Ivan der Marel fc Franxl 
1 19931 : Ide Lorenzi et al' 20081 ) and the luminosity itself (cor- 
responding to ho) within the field-of-view, for a total of 
Nkin ~ 365 kinematic observables (73 Voronoi bins times 
5 sets of Gauss-Hermite coefficients). These are used to con- 
strain the particle system as a whole, without distinguishing 
bf between NCD and NCS populations. The particle model 
is seeing convolved with a Gaussian PSF having a FWHM 
of 0'.'23 ( 5 pc) by means o f the Monte-Carlo method pre- 
sented in Ide Lorenzi et al.l (2008): When the model kine- 
matics are computed, each particle is temporarily replaced 
by Npp pseudo particles with randomly selected positions 
having probabilities given by the PSF. In the present work, 
we adopted Npp — 5. 



3.2 Initial conditions 

We set up spherical initial conditions using the major- 
axis density profile of the NCS. Setting M/Lk ~ 1, the 
mass density is normalised to unit mass and the self- 
consistent grav itational potential is computed. Following 
I GerhardI (|l99ll ). the isotropic distribution function is com- 
pu ted and used to generate a se t of equal mass particles as 
in iDebattista fc SellwoodI (|2000l '). Finally, 30% of the par- 
ticles are randomly assigned to the "disc" population. We 
produced two realisations of these initial conditions, a low 
resolution one with 0.75M particles which allowed us to ex- 
plore parameter space quickly and a higher resolution ver- 
sion with 6M particles. 



3.3 Gravitational potential 

The M2M method works by adjusting the weights of indi- 
vidual particles while they are evolved along their orbits. 
For this orbit integration the gravitational potential of the 
system is needed. This section details the methods used to 
construct the dynamical models presented in Section |4] be- 
low. In brief, low resolution models (Ml) are built using an 
FFT method for the potential. Then high resolution mod- 
els without (M2) and with intermediate mass black holes 
(IMBHs) (M3-M8) are built with the pot ential compute d 
on a spherical mesh of spherical harmonics (|Sellwoodll2003l ) . 

We assume that the mass distributions of the MD, NCD 
and NCS follow their luminosity distributions Ji (i G {MD, 



NCD, NCS}). For a constant mass-to-light ratio M/ Li, the 
corresponding mass density is pi — {M/Li)Ji, with Ji the 
deprojection given in Section r2.1.1l 

The potential in models Ml was obtained on a Carte- 
sian grid. Each Cartesian grid consists of TV = 128^ grid 
cells. The NC grid extends to ±50 pc along each direction, 
whereas the MD grid extends to ±3000 pc. This allows us 
to resolve the scale-heights of the MD and the NCD and the 
half- mass radius of the NCS. The potential on each grid is 
calculated using the Fourier convolution theorem. We assign 
to each mesh point a mass from the corresponding density 
distribution. The potential is then obtained by a convolu- 
ti on with the Greens function. We employ the FFT method 
of IPress et al.l 11992) to perform the convolution. 

We pre-compute the individual gravitational potentials 
<I>f^^ generated by Ji, for unit mass-to-light. This proce- 
dure allows us to quickly obtain the total gravitational po- 
tential for any choice of M/Li through a weighted sum 
<1> = M / Li^f'^ , which is kept constant for each model 
run. 

For a modelling run, we initially tabulate "5^ = 
M/Li^f'^ on individual Cartesian grids. Forces at grid 
points are computed by finite differences. Individual par- 
ticle acc elerations are then approxim ated by a cloud-in-cell 
scheme (|Hocknev fc Eastwood|[T988l ) to interpolate the grid 
point forces to the particle position. 

For models in the series M2-M8 w e use a spherica l 
harmonic potential solver as described in ISellwoodI l|2003h . 
We use the disc and spheroid particle populations (particle 
weights are converted to mass via the associated M/Lk) 
to calculate the gravitational potential of the entire parti- 
cle system. Thus, the contribution of the MD is neglected. 
As discussed at the end of Section 12.1.11 the error associ- 
ated with this approximation is expected to be very small. 
For models using the spherical harmonic potential solver, 
the potential is updated after every M2M correction step 
(and temporally smoothed) . We use potential expansion co- 
efficients up to Imax = 8 with 300 (unevenly spaced) radial 
bins to Vmax = 200 pc. The width of the innermost bin 
is 0.2 pc and of the outermost bin is 4.2 pc, which is still 
smaller than the FWHM of the PSF. 



3.4 Resampling a particle model 

We use the final particle dynamical model as initial con- 
ditions for Ai'-body simulations. In order to do this, it is 
best that the particles have a narrow range of masses; this 
ensures both higher effective mass resolution and a lower ar- 
tificial two-body relajcation rate. The models are t herefore 
built using the re-sampling technique described bv lDehn"enl 
(2009). This section closely follows Dehnen (200^ work. We 
generate the models using a fiat weight prior w — . 
The particle models are re-sampled every 100 M2M correc- 
tion steps if the ratio of largest to smallest particle weight is 
max{i(Ji}/ minlwi} > 10. Because we do not normalise total 
weight (it is only constrained by the observables, in partic- 
ular by yloo), the weight of a re-sampled particle is set to 
uifc = '^f,Wi^oid- The phase-space coordinates (xfc,Vfe) 
of the k*** re-sampled orbit are set to the i*" ori ginal trajec- 
tory if 



a < 7(fc - 1/2) < a+u i,ke [1, N] 



(4) 
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with mean relative normalised weight 7 = N^^ ^^"fi, cu- 
mulative relative normalised weight d — '}2k<i 'Y* ^^'-^ 
7i = Wi/w the relative weight. Orbits with 7^ < 7 are re- 
sampled at most once, whereas orbits with 7^ > 7 produce at 
least one copy. If a trajectory gets re-sampled, the first copy 
gets the phase-space position of the original particle. For ad- 
ditional copies we randomise positions (x±d, j/±d, z±d) with 
d = (x^ + j/^)^^^/100 exploring the eight distinct combina- 
tions of plus and minus signs. The position of any additional 



copy is set to (a; -I- r^d, y + ryd, z + r^d) with Tx 



'1,1] 



uniform random numbers. We do not alter velocities and ev- 
ery copy keeps the velocity of the original particle. At small 
radii, r < 0.2 pc, we rotate the particles randomly about the 
z-ajds otherwise the resulting orbits are too closely spaced. 
Our implementation of resampling conserves total particle 
number but not individually for the NCD and NCS. When a 
particle is re-sampled more than once all daughter particles 
inherit its affiliation to either the NCD or the NCS. 

Particle re-sampling at work is illustrated in Figure [2] 
The top panel compares a final particle weight distribution 
of a model generated using nmagic without re-sampling 
with the peaked distribution of a corresponding model built 
including re-sampling. The model with a narrow weight dis- 
tribution has a higher effective resolution and hence suf- 
fers less from shot noise than its counterpart with a broad 
distribution. Starting from a spherical particle population 
the NCD shown in the bottom panel is obtained by com- 
bining NMAGIC with re-sampling. The eff ective number of 
parti cles is defined as N^ff/N = uP' j-uP- (Ide Lorenzi et all 
[2003) • Then, the particle models shown in Figure [5] have 
1 - N^ff/N = 0.67 without re-sampling and 1 - N^ff/N = 
9.8 X 10~^ with re-sampling, i.e. re-sampling leads to ~ 3x 
higher N^ff. 



3.5 Run parameters 

Based on our experience in previous work, we set the value of 
the force-of-change parameter e = 8 x 10~* and the temporal 
smoothing parameter a = 2.1e. We set the entropy param- 
eter ^ = 2 x 10~ ^ . (The par ameters (e, a, /i) a r e defi ned in 
ISver fc Tremaind (| 19961 ') and Ide Lorenzi et all l|2009l ).'l We 
used timesteps St = 673 years, with x^-M2M correction 
steps every 20 timesteps. For comparison, a circular orbit 
at r = 13.4 pc, which contains half the particle mass, takes 
1.82 Myr. 



4 DYNAMICAL MODELS 

In this section we construct dynamical models for the cen- 
tral region of NGC 4244 to assess its intrinsic kinematics 
and to constrain the NC mass. We investigate axisymmetric, 
two-component models for different combinations of mass- 
to- light ratios, fitting the photometry and NIFS integral field 
kinematic data. 

All models are constructed including re-sampling of par- 
ticle coordinates as described in Section [3.41 In order to re- 
duce computational cost, the bulk of the modelling is per- 
formed using 0.75M particles, with the gravitational poten- 
tial calculated on the Cartesian grids and held fixed through- 
out. In these models trajectories are integrated with a stan- 
dard leapfrog scheme with a fixed time step. The models 



o 



a- O 
2 O 
O 



O 
O 



-1 -0.5 0.5 

Log(w/wQ) 



u 

^ 



-10 - 



o 

■3 



-10 ■ 




Figure 2. Top: comparison of the particle wciglit distribution 
of model M2 generated including re-sampling (solid line) and a 
corresponding model produced without re-sampling (dot-dashed 
line). Bottom: Distribution of a subset of NCD particles of the 
M2 model, generated using NMAGIC with re-sampling. The three 
panels show different projections of the particle coordinates. Note 
that the initial NCD model was spherical. 



are constructed in a two step process. First, we start with 
the spherical isotropic 0.75M particle model and evolve it 
using NMAGIC to generate a particle realisation with de- 
sired luminosity distribution, fitting simultaneously but sep- 
arately the NCS and NCD photometric constraints. Since 
only photometric (not kinematic) constraints are fit at this 
stage, the velocity scaling is arbitrary and only the ratio of 
Af/LNCD and M/Lncs matters for the shape of the gravita- 
tional potential. In order to compute the gravitational po- 
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tential the ratio (M/Lncd)/(M/Lncs) is fixed at 0.2/1.8. 
The resulting model then serves as a starting point to si- 
multaneously fit both the photometric and kinematic con- 
straints {Noba ~ 2365 observables) for diflerent combina- 
tions of M I I/NCD and M / -Lncs ■ During this adjustment 
phase, we typically evolve for 4M timesteps (2.7 Gyr) and 
apply 200K M2M correction steps. The particle system is 
then relaxed for a further lOOK timesteps (67 Myrs) with- 
out changing particle weights. We refer to the final models 
as series Ml. We vary A//Lncs between 1.0 and 4.8, whereas 
we use values for M/Lncd of 0.1, 0.2 and 0.4. The influence 
of the mass-to-light ratio of the MD on the quality of the 
fit is expected to be negligible; we therefore keep it constant 
at M/Lmd = 0.7 compatible with estimates from integrated 
colours. 

The results are presented in Figure [S] which illustrates 
how the quality of the model fit changes with mass-to- 
light ratios. The top panel shows Ax^ = ~ minfx'^}) 
whereas AxLn of the kinematic observables alone (Ax'^ 
"marginalised" over the ylj^'s) is shown in the bottom panel. 
The 68% confidence limits (la) are computed as \/2Nobs 
following Ivan den Bosch fc van de Ve n (2009). These limits 
are Ax^ = 68.8 and Axtin = 27.0, respectively. 

Regardless of whether AxLn or Ax^ is used, the result- 
ing best model has mass-to-light ratios M/Lncd = 0.4 and 
M/Lncs = 1.4. Note that the models with M/Lncd = 0.1 
or M/Lncd ~ 0.2 and M/Lncs ~ 1.6 reproduce the 
NIFS data with comparable quality and are in better agree- 
ment with M/Lncd = 0.1 — 0.25 estimated from integrated 
colours (see Section [2.ip than is M/Lncd ~ 0.4. Moreover, 
M/Lncd =0.1 als o agrees with t he estimates from Bruzual 
& Chariot models (|Bruzual fc Chariot 200a) . 

The range of acceptable NCS masses is estimated as 
the range for which Ax^ (AxLu) is below la confidence 
after "marginalising" over M/Lncd. Although Ax^ is the 
more appropriate quantity to discriminate between models 
because both photometric and kinematic constraints are im- 
posed on the models, the allowed NCS mass ranges deter- 
mined using either Ax^ or AxLn agree with each other. We 
obtain a NCS mass of Mncs = l-6to.2 10^ within 
~ 42.4 pc. This is almost an order of magnitude higher 
than the lower limit of ~2.5xl0^ M0 obtained from the 
observed velocity of an Hll re gion at a project ed distance 
of 19 pc from the NC center (|Seth et al.l 12006*). The mass 
within ~ 15 pc is 1.0 x 10^ M©, which agrees with the mass 
within the same radius obtained from the JAM models in 
Hll, (1.1 ± 0.2) X lO'^ Mq. The mass of the NCD is not as 
well constrained: we obtain 3.6 x lO'^ Mg for M/Lncd = 0.1 
and 14.4 x 10^ Mq for M/Lncd = 0.4. 

At this point models with 6M particles are constructed 
starting from the spherical isotropic 6M particle initial con- 
ditions. We start by generating high resolution versions of 
models Ml with A//Lncd = 0.1. We again use photomet- 
ric followed by photometric+kinematic constraints. But now 
we replace the FFT method with a spherical harmonics po- 
tential solver in order to obtain higher spatial resolution 
at the centre. These models also use a Runge-Kutta time- 
integrator with ad aptive timestep (using routine ODEINT of 
iPress et al.l (|l992l ) with accuracy parameter 10 ®) to allow 
a comparison with the models which include a black hole, 
presented in the next section below. Since the potential is 
computed via an expansion in Z, m spherical harmonics anal- 
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Figure 3. Ax^ = ~ rnin{x^} as a function of NCS and 
NCD /f-band mass-to-light ratio. Upper panel: Total Ax^ for 
the model fit to the photometric and kinematic target observ- 
ables. Bottom panel: ^Xkin °^ kinematic observables only. 
Dashed horizontal lines corre s pond to 68% (Icr) confidence after 
Ivan den Bosch &: van de Venl l l2009l) . Symbols are as follows: The 
Ml series for M/Lncd = 0.1, 0.2 and 0.4 are shown as open 
(red) circles, (black) squares and (blue) triangles, joined by dot- 
dashed, solid and dotted lines, respectively. Models in series M2 
are indicated by the solid circles joined by the red dot-dashed 
lines. 



ogous to the photometric constraints, we include the same 
terms in the expansion of the potential as for the luminosity 
density {i.e. non-zero Af^^ and AfJ^^, cf. also Section [3. 3p . 
The potential is recalculated after every M2M correction 
step. We refer to the resulting models as series M2. 

The 6M particle models illustrate several interesting 
points. The NCS mass estimated using 6M particles agrees 
with the estimates presented above using the Ml models. 
This suggests that the inferred NCS mass is robust with 
respect to how the models have been constructed, in partic- 
ular to the potential solver, integration scheme and number 
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of particles. Increasing the number of particles from Ml to 
M2 decreases Xfcin of the corresponding best models by only 
a small amount (if at all) with respect to the Axtin confi- 
dence limit. This indicates that the model fit to the NIFS 
data is dominated by the uncertainties in the data while the 
contribution of shot noise to Xkin is negligible. On the other 
hand, is reduced considerably mainly due to a reduction 
in Xa;™EI- This suggests that XAim is dominated by Pois- 
son noise. Generally XAim < NAim because of the temporal 
smoothing, and both the Ml and M2 models reproduce the 
photometric data very well. 

We present the intrinsic kinematics of the best-fit model 
in series M2 in Figure |4] We computed the radial, tangential 
and vertical dispersions, (Ju, ctv and respectively, and plot 
the anisotropics f^^p — 1 ~ (o-„/cr„)^ and /3z = 1 — (cru,/(Tu)^. 
In agreement with Hll we find that /3z < 0. Also Vrms has a 
central minimum. Hll found that the combination of these 
two properties provide important constraints on the amount 
of mass that the nuclear cluster could have accreted in the 
form of star clusters, as we shall also see below. 



4.1 Adding intermediate mass black holes 

Some galaxies a re known to harbour both a NC and a mas- 
sive black hole (ISeth et al.ll2008al : iGraham fc Spitled [2OO9I : 
iNeumaver fc Walcherll20li l. Since mass and anisotrop y are 
degenerate with each other l|Binnev fc MamonI 1 19821 ') , we 
wish to explore how adding an IMBH might change (3^. We 
therefore also generate models including an IMBH at their 
centre with the aim of finding a robust upper IMBH mass 
limit compatible with the observations. We construct these 
models in a manner analogous to models M2 above, using 
a spherical harmonics potential solver in order to obtain 
higher spatial resolution at the centre. These models again 
use a Runge-Kutta time-integrator with adaptive timestep 
for higher accuracy i n the vicinity of t he IMBH. Using two- 
integral JAM models l|CappellariH2008l ') Hll obtained an up- 
per limit of M. < 10^ M© on any black hole that may 
be present. We revisit this estimate with our more general 
three-integral modelling. 

The IMBH is represented by a Plummer potential with 
scale-length set to 0.02 pc. We generate models M3 to MS 
for various IMBH mass fractions. The IMBH mass fractions 
C = M./{M/Lncs X IO^Mq) are given in Table (2] For each 
series of models the black hole mass is then given as M. = 
( X M/I/NCS- The results of the models are presented in 
Figure O which illustrates how the quality of the model fit 
changes with A//Lncs and The left panel shows Ax^ = 
X^ ~ min{x^}, whereas AxLn of the kinematic observables 
alone is shown in the middle and right panels. The 68% 
confidence limits are as in the above i.e. Ax^ = 68.8 and 



AXfci. 



27.0. 



As expected, the minimum x^ along a given line in Fig- 
ure [5] shifts towards smaller AZ/Lncs with increasing IMBH 
mass fraction (see especially the right panel of Figure [S]). 
Each line in Figure [5] intersects the confidence limit (dashed 
horizontal line) twice (in the case of series M3 and M4 the 
modelling sets need to be extrapolated). The intersection 
with larger M/L-^cs corresponds to the largest admissible 



^ We use the same MC Ai^ errors as for models Ml. 




0.2 0.3 

R (arcsec) 

Figure 4. The intrinsic kinematics of the NC in the best-fit model 
in the M2 series. From bottom up we show the second moment 
of line-of-sight velocity Vrms, line-of-sight velocity V, the circu- 
lar velocity Vc, the vertical anisotropy Pz, and the tangential 
anisotropy /30. The top row shows the Toomre-Q of the NCD 
only. 



IMBH mass along each line. The largest IMBH mass com- 
patible with the data would be obtained at a line that inter- 
sects the horizontal line only once, at its minimum x^- Using 
Axiin shown in the right panel of Figure (5] the dashed line 
with A^/Lncs = 1.2 leads to the IMBH mass upper limit 
of M. = 5.0 X 10^ Mq. This upper limit is larger than the 
one found in Hll using JAM models, presumably reflecting 
the greater orbital freedom presented by 3-integral versus 
2-integral models. 

If we use Ax^ shown in the left panel of Figure[S]instead 
of Axiin, w6 find that even larger IMBH masses (up to a fac- 
tor 2 or higher) are compatible with the data. Nonetheless, 
we use the more conservative IMBH mass range provided by 
the NIFS data alone. 

Figure [6] shows a comparison of the best fit models in 
series M2 and M6 (indicated by the orange stars in the 
middle and right panels of Figure [S} with the integral field 
NIFS kinematic. The model fits to the NIFS data are ex- 
cellent. Xkin for the best models M2 and M6 are 194.715 
and 207.532, respectively. Figure [7] shows a comparison of 
the best fit models in the M2 and M6 series with major- 
axis kinematic data extracted from the NIFS data. The 
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Figure 5. Ax^ = — niui{x^} as a function of NCS mass-to-ligiit ratio for various IMBH mass fractions Q = yL,/{M/ LffQg). Left panel: 
Total Ax^ for ttie model fit to the photometric and kinematic target observables. Middle panel: ^Xkin °^ kinematic observables 
only. Right panel: AxLv. as a func tion of M, instead of M/Lncs- Dashed horizontal lines correspond to 68% (Icr) confidence after 
Ivan den Bosch Sc van de VenI l l2009t) . In the left and middle panels the symbols are as follows: The M2 series is shown as solid (red) 
circles joined by red dot-dashed lines. Models including a central IMBH are shown as open squares joined by (black) solid (M3), (green) 
dashed (M4), (blue) dot-dashed (M5) or (orange) dotted (M6) lines. The open triangles are for series M7 and M8. The M7 and M8 
series are joined by (pink) solid and (grey) dashed lines respectively. IMBH mass fractions for models in the MS to M8 series are given 
in Table [2] In the right panel the different symbols represent: circles joined by (black) dot-dashed line: M/Lncs = li squares joined 
by (red) solid line: M/Lncs = l-li squares joined by (green) dashed line: AZ/Lncs = 1-2, squares joined by (blue) dot-dashed line: 
A^/^NCS = squares joined by (pink) dotted line: M/L^jcs = 1-4, triangles joined by (grey) solid line: M/Lncs = (single) crossed 
circle: M/Lncs = 1-6, (single) dotted circle: M/Lncs = 1-7, (single) cross: M/Lncs = l-9i and (single) diamond: M/Lncs = 2.2. The 
open and solid (orange) stars in the middle and right panels mark the best fit M2 and M6 models, respectively. 



Series 




M/Lncd 


M/Lncs 


^ (M/Lncs )x 105 Mq 


Pot. 


Int. 


Ml 


0.75 


0.1,0.2,0.4 


1.0 - 4.8 




FFT 


Leapfrog 


M2 


6 


0.1 


1.0 - 2.2 




Sph. harm. 


RUNGE-KUTTA 


M3 


6 


0.1 


1.1 - 1.5 


0.6 


Sph. harm. 


Runge-Kutta 


M4 


6 


0.1 


1.1 - 1.5 


0.9 


Sph. harm. 


RUNGE-KUTTA 


M5 


6 


0.1 


1.1 - 1.5 


1.2 


Sph. harm. 


RUNGE-KUTTA 


M6 


6 


0.1 


1.2 - 1.4 


2.3 


Sph. harm. 


RUNGE-KUTTA 


M7 


6 


0.1 


1.0 - 1.3 


3.5 


Sph. harm. 


RUNGE-KUTTA 


M8 


6 


0.1 


1.0 - 1.4 


5.8 


Sph. harm. 


RUNGE-KUTTA 



Table 2. Summary of the modelling runs. Columns from left to right are: Name of the series of models, number of particles, mass-to- 
light ratio of the NCD, mass-to-light ratio of the NCS, IMBH mass fraction (J, the method used to compute the potential and the time 
integration scheme. 



model kinematics are computed from a Gauss-Hermite fit 
to the line-of-sight velocity distribution in the corresponding 
Voronoi bins. To compute the temporally smoothed LOSVD, 
we use 27 bins in velocity within a range of width 300 
km s~^, centred on the corresponding NIFS line-of-sight ve- 
locity. Note that we did not fit the full LOSVD itself, in- 
stead we constrained the particle models using luminosity- 
weighted moments as described in Section [3. II 



5 AT-BODY SIMULATIONS 

Our best dynamical model without an IMBH is M2 with 
M/Lncs ~ 1-5. After building this model we used it for a 
number of A^'-body experiments testing its stability and evo- 
lution by using it as initia l conditions. The simulations, were 
evolved with pkdgrav l|Stadell 1200 ll ). an efficient, multi- 
stepping, parallel treecode. In all cases we use an open- 
ing angle 9 = 0.7. We used base timestep At — 0.1 Myr 



and changed timesteps of individual particles such that 
St = At/2" < rj{e/ay^^, where e is the softening and a is the 
acceleration of the particle, with n as large as 29 allowed. 
We set rj = 0.03, a quite conservative value. 



5.1 Stability test 

In constructing our dynamical model we have assumed that 
the NC in NGC 4244 is very likely axisymmetric. Hll found 
no evidence of non-axisymmetry in the NC of M33. Its PA 
is consistent with that of its main disc and its apparent 
ellipticity is consistent with a vertical flattening of g = 0.7, 
the average observed in the NCs of edge-on late- type galaxies 
(Seth et al. 200^). There is also only a small misalignment 
between photometric and kinematics major axis. At present, 
M33 is the only galaxy in which the axial symmetry of the 
NC can be determined. 

Our first A'^-body simulation therefore tests the stability 
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Figure 6. Symmetrised NIFS integral-field kinematic data within ±0.7" of the central region of NGC 4244 (top row) compared with 
corresponding luminosity- weighted data extracted from the best-fit models in series M2 (middle row) and M6 with M/Lj^qj^ = 0.1, 
Af/LjMcg = 1.3 and Mbh = 3.0 X 10^ Mq (bottom row). These models are indicated by the open and solid orange stars, respectively, 
in the right panel of Figure [5] From left to right arc shown: line-of-sight velocity v, line-of-sight velocity dispersion a and higher order 
Gauss-Hermite moments h3 and h4. 



of the model against non-axisymmetric perturbations, par- 
ticularly the bar instability, which plagues rapidly rotating 
systems. After evolving the best fit model in series M2 for 50 
Myr (the rotation period at 5 pc being 0.63 Myr), the model 
remained axisymmetric with no hint of a bar or spirals. The 
top panel of Figure |4] plots the Toomre-Q = (7uk/(3.36GE) 
of the NCD, where k is the epicyclic radial frequency, G is 
the gravitional constant and E is the surface density. The 
stability of the system stems from the high Toomre-Q of the 
NCD, which is everywhere greater than 10. 



5.2 Accretion simulations 

Hll explored the hypothesis that NC formation is a result of 
star cluster (SC) accretion. They showed that the observed 
kinematics of the NC in NGC 4244 are not consistent with 
accretion of more than ~ 50% of its mass in the form of 
SCs. Specifically, once the accreted mass fraction exceeded 
this value the resulting central Vrma was no longer a mini- 
mum. We will see here too that Vrms starts to lose its central 



minimum once the accreted mass fraction becomes too large, 
and our constraint is even more stringent than that of Hll. 
On the other hand, Hll showed that the negative /?z of the 
NC was only possible if it accreted > 10% of its mass as 
SCs on highly inclined orbits. Those simulations assumed 
SCs accreting onto a pre-existing NC, either as a NCD or as 
an isotropic NCS. Therefore we next subject our best-fitting 
model of the NC to SC accretions. 

In order to permit the model star clusters to sink to 
the centre via dynamical friction, we introduce the best fit 
model in the M2 series inside a particle main disc with an 
exponential profile, because NGC 4244 is a late-type, bulge- 
less galaxy. We use the same model for the MD as did Hll, 
I.e. four million multi-mass particles with masses ranging 
from 7M0 within the inner 20 pc increasing to 1.2 x lO'^ Mq 
in the disc's outskirts. The distributions of masses and soft- 
enings of the main disc particles are shown in Figure 6 of 

1/3 

Hll; the softening is related to particle mass via ep oc irip . 

We accrete three of the model SCs described in Hll 
which we term Gl, G2 and G3 in order of increasing mass. 
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Figure 7. Comparison of the best-fit models in series M2 and 
M6 (indicated by the open and solid orange stars in Figure (SJ 
to the kinematic data along the major-axis extracted from the 
NIFS data. Models and data are shown as hnes and open circles, 
respectively, with M2 indicated by the dot-dashed (red) line and 
M6 by the dotted (black) line. From bottom to top are shown: 
velocity v, velocity dispersion a, and Gauss-Hermite moments /i3 
and hi. The model kinematic data were computed using a Gauss- 
Hermite fit to the velocity distribution along the line-of-sight. 



Their properties are listed in Table [S] As in Hll, we define 
the concentration of each star cluster as c = \og{Ras/Rc) 
where ReS is the half mass radius (effective radius) and R^ 
is the core radius, where the surface density drops to half 
of the central. These were allowed to accrete onto the NC 
starting from circular orbits at 50 pc from the centre. We 
start the SCs from 4 different inclinations relative to the 
NCD: 0°, 30°, 60°, and 90°. The SCs require about 40 Myrs 
to accrete onto the NC. 

The results of these accretion simulations are shown in 
Figure[Hl The low density star cluster Gl is disrupted at ~ 8 
pc from the centre. Thus it barely perturbs the kinematics 
of the NC. However we note that there is a general tendency 
for pz to increase slightly within RcS, suggesting that even 
this mild ~ 1% mass accretion can alter the kinematics. The 
more massive, denser, clusters G2 and G3, both of which sink 
all the way to the centre, perturb the kinematics much more. 
In all cases Pz increases; for G2 averaged within R^s is 
nearly zero, while /3z > everywhere within iioff when G3 
is accreted. 

Accreting G3 (which has a mass ~ 13% that of the NC) 
also raises the central Vrms. Although Vrms at the centre 
remains smaller than at ReS, there does not seem to be much 
room for further significant accretions without making Vrms 
centrally peaked, unlike the observations. This is likely to 
hold also of this mass fraction arrives as many smaller star 
clusters, provided that the star clusters are dense enough 
that some fraction of their stars survives all the way to the 
centre of the NC. 



Model 


M, 


^off 


c 


Hll name 




[xIO^Mq] 


[pc] 






Gl 


2 


1.11 


0.12 


C4 


G2 


6 


1.11 


0.16 


C5 


G3 


20 


2.18 


0.12 


C3 



Table 3. The star clusters used in the accretion simulations. M» 
is the stellar mass of the SC, iJoff is the effective (half-mass) 
radius, and c is the concentration (defined in the text). For com- 
parison, the last column lists the name used for the model in 
Hll. 




1000 r n 
10^ 10'° 10" 10^^ 

"goI (Mo) 

Figure 9. Mass of the NC of NGC 4244 versus the mass of 
the host galaxy comp ared with the Mqmo ~ ^gal relations of 
iFerrarese et al.l ll2006h . The solid red and black lines show the cor- 
relations for NCs and SMBHs in early-type galaxies, respectively, 
with 1 a confidence levels shown as dashed lines. The blue lines 
show the relation of lScott &: GrahamI l l2012h . 

6 DISCUSSION AND CONCLUSIONS 

We have performed a dynamical study of the nuclear cluster 
in the edge-on spiral galaxy NGC 4244 taking into account 
different morphological components, which are the galaxy 
main disc (MD), nuclear cluster disc (NCD) and nuclear 
cluster spheroid (NCS). We have constructed axisymmetric 
dynamical particle models accounting for the MD, the NCD 
and the NCS. We find a total NCS mass Mncs = x 
10' ^ Mp) with ir i appro ximately 42.4 pc ( 2"). B oth the fits 
of ISeth et all (12005a') and of iFrv et all (|l999l ) show that 
there is no obvious bulge com ponent in NGC 424 4. Using 
the 2MASS'L&xge Galaxy Atlas l|jarrett et al.ll2003l ) JT-band 
magnitude, the total luminosity of the galaxy is 3.2 x 10^ Lq, 
and thus the galaxy stellar mass is 2 x 10® M©. Figure [5] 
plots the NC mass compared with the Mcmo — Mgai rela- 
tion. The NC sits above this relation. 



6.1 Vertical anisotropy 

The kinematics are moderately tangentially anisotropic in- 
side Reft with an anisotropy parameter Pz ~ —0.1. This is 
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Figure 8. The effect of accreting star clusters onto the NC in NGC 4244. From left to right these show the effect of accreting SCs Gl, 
G2 and G3. The black, red, green and blue solid lines show SCs accreting from 0°, 30°, 60° and 90°, respectively. The dashed line shows 
the initial NC. 



in good agreement with the JAM models of Hll. Hll showed 
that /3z < requires high-inclination infall of star clusters 
onto a pre-existing nuclear cluster. In our accretion simula- 
tions onto a more realistic model of the NC we found that 
even the accretion of a star cluster of just 13% the mass is 
enough to erase the vertical anisotropy. This raises questions 
about whether such anisotropy can be due to accretion at 
all. It also hints that, unless we are observing the NC of 
NGC 4244 at a special time, it cannot sustain accretion of 
> 10% mass as suggested by Hll. 

We therefore tested whether the assumption of a per- 
fectly edge-on nuclear cluster may bias the modelled ver- 
tical anisotropy to negative values if the real inclination 
is somewhat smaller. The smallest inclination at which we 
were able to deproject the NCD photometry was 83°. Us- 
ing a model deprojected at this assumed inclination and 
the observed kinematics, we built nmagic models assum- 
ing M/Lt^cs = 1-5 starting from the best-fit edge-on model 
M2. The dashed black lines in Figure [TO] show that the 2- 
D anisotropy, — 1 — a^/aj^, and 3-D anisotropy, Bz = 
1 — 2(j1/{a\ + f^), are barely changed compared to the 
edge-on case (solid black lines) and remain negative. Thus 
a negative vertical anisotropy is not an artifact of assuming 
the nuclear cluster is perfectly edge-on. 

We finally explore whether the recovered /3z changes if 
we include IMBHs in the models. In Figure [TOl we plot both 
the vertical anisotropies for vary M, . While increasing M, 
raises the vertical anisotropy, it still remains negative within 
ReS- We conclude that the NC must be vertically anisotropic 
even if a black hole were present. 



6.2 Summary 

We have built dynamical models of the NC in the nearby, 
edge-on late-type galaxy NGC 4244. Using particle re- 
sampling, we were able to obtain a narrow distribution of 
weights in our NMAGIC models allowing us to use the mod- 
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Figure 10. Profiles of 2-D (top) and 3-D (bottom) vertical 
anisotropy. The black, blue, green and red lines correspond to 
M. = 0,0.8,1.23, and 3.0 X lO^^ M©, respectively The dashed 
lines correspond to assuming that the nuclear cluster is inclined 
at 83° instead of being perfectly edge-on. 
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els as initial conditions in A'^-body simulations. Our results 
can be summarised as follows: 

• We find a mass of the spheroidal component of the 
NC, Mncs = 1-6±S5:2 X 10^ Mq within 42.4 pc. The mass 
within 15 pc is ~ 1.0 x IO^ Mq, in very good agr eement 
with the value estimated by iHartmann et al.l (|201ll ) using 
two-integral JAM models. This mass puts the nuclear cluster 
above the MNC-Mcai relation. 

• The mass of the bluer disc component of the nuclear 
cluster is less well constrained and covers the range 3.6 x 
lO^Mo < Mncd < 14.4 X IO^Mq. 

• Our three-integral models are consistent with no black 
hole as well as with a black hole at least as massive as 4.6 x 
10^ Mq . This upper limit is larger than the one allowed by 
two-integral JAM models. 

• Simulations show that the model without a black hole 
is stable against axisymmetric perturbations. This stability 
derives from the large Toomre-Q of the system. 

• Regardless of whether a black hole is present or not, 
red and of whether the nuclear cluster is perfectly edge-on 
or not, and are both negative. Accretion of a star 
cluster of as little as 13% by mass is enough to drive /3z to 
positive values, regardless of the orbital geometry. It remains 
unclear, therefore, how /?z < arose. 
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